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We study gravitational waves from a particle moving around a system of a point mass with a 
disk in Newtonian gravitational theory. A particle motion in this system can be chaotic when the 
gravitational contribution from a surface density of a disk is comparable with that from a point mass. 
In such an orbit, we sometimes find that there appears a phase in which particle motion becomes 
nearly regular (so-called "stagnant motion" or "stickiness" ) for a finite time interval between more 
J"""""- \ strongly chaotic phases. To study how these different chaotic behaviours affect observation of 

gravitational waves, we investigate a correlation of the particle motion and the waves. We find that 
■ such a difference in chaotic motions reflects on the wave forms and energy spectra. The character of 

(N ' the waves in the stagnant motion is quite different from that either in a regular motion or in a more 

^q, strongly chaotic motion. This suggests that we may make a distinction between different chaotic 

r—t • behaviours of the orbit via the gravitational waves. 
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INTRODUCTION 

cr. 

Chaos appears universally in nature and it is expected to be a fundamental tool to understand various nonlinear 
phenomena. Historically, research on chaos started from the famous three-body problem by Poincare, and so far many 
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attempts to reveal the character of chaos have been done in dynamical systems. Following such works, considerable 

research in Newtonian gravity and general relativity have been done 0, [nl q q, q h, 

Q, 0, Q Q, H 

\2l\ . But most of it, especially work on relativistic systems, has discussed only whether or not 



. chaos occurs mainly by using the Poincare map and the Lyapunov exponent. However, we know there appear various 



types of chaotic behaviours depending on the strength of chaos and its analysis will play a very important role to 

o 
t> 
o 



' understand the essence of nonlinear dynamics [lcl . A coquet view of how one can extract and use information 



from a chaotic system may be also missing. So one can address two new important issues in the research of chaos in 
• , Newtonian gravity and general relativity One is, of course, to make clear the character of chaos for each system, and 
■ the other is to find some methods to extract useful information from chaotic systems. As for the first point, we have 



recently shown the possibility to classify the character of chaos in a system of a spinning particle moving around a 
Schwarzschild black hole 2l|. The method used in 2l| is a power spectrum analysis of the particle orbit. The spectra 



are mainly classified into a power-law type and a white-nois e type. As a result, we find that there is a close relation 
between the so-called "stagnant motion" (or "stickiness" \ut lisl . I20I]) and a "power-law" spectrum. 

As for the second point, an indirect method to extract information from chaotic systems is required for the following 
reason: In chaotic systems in astrophysics, it is sometimes too difficult to observe chaotic motions directly. Because 
these systems are often far from the earth and the ambient surroundings of these systems may not be clean. Therefore, 
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in this paper we propose the use of gravitational waves as a new method to analyze chaos. The reason we choose 
gravitational waves is as follows: Detection of gravitational waves is one of the greatest challenges in experimental 
and theoretical physics in this decade. Several kilometer-size laser interferometers, such as TAMA 22j, LIGO 23], 



and GEO [24| are now in operation. In addition to these ground-based detectors, the Laser Interferometer Space 
Antenna (LISA) with an arm length of 5 x 10 6 km has been proposed and is planned to start observation in the near 
future [251 ]. Gravitational waves will bring us various new information about relativistic astrophysical objects. If 
we detect gravitational waves and compare them with theoretical templates, we may be able to determine a variety 
of astrophysical parameters of the sources such as their direction, distance, masses, spins, and so on. The direct 
observation of gravitational waves could resolve strong-gravitational phenomena such as a black hole formation. 
Furthermore, we may be able not only to verify the theory of gravity but also to find new information at high density 
or to recover new physics in a high energy region. 

In [l7|, we analyzed the gravitational waves from a spinning test particle in a Kerr black hole. We find that there 
is a difference between the spectra of the gravitational waves from a chaotic orbit and from a regular one. There 
appear many small spikes in the spectrum of the chaotic orbit. However, as we mentioned, there are various types of 
chaotic motions, and it is important in the analysis of such a dynamical system to know which type of chaos appears 
as well as to show the difference between a regular motion and a chaotic one. Hence, in order to study whether one 
can make a distinction between various types of chaos by use of gravitational waves, we should reanalyze them in a 
chaotic system. 



26, 



27|. 



As a concrete model of a chaotic system, here we consider a point mass with a thick disk in Newtonian gravity 
This model mimics a system of a black hole with a massive accretion disk 2^] . Saa analyzed this system and showed 



that a particle motion is chaotic 



271 ]. This model can describe almost regular to highly chaotic motion by changing 



the ratio of a disk mass to a black hole mass. In 27], however, only the Poincare map has been analyzed to judge 
whether chaos occurs or not, and the characteristics of chaos have not been studied much. 

So our strategy in this paper is the following: First, we analyze the particle motion and make clear the characteristics 
of chaos appearing in this system. Secondly, we evaluate the gravitational waves from such a system by use of the 
quadrupole formula. Finally, to study some observational feature of chaos appearing in the gravitational waves, we 
investigate correlation between types of chaotic motions and gravitational waves, and then point out a possibility to 
extract information about this chaotic system from the gravitational waves. 

This paper is organized as follows. In Sec. [TTJ we shall briefly summarize the basic equations. Numerical analysis 
results will be presented in Sec. IIIIl The summary and discussion follow in Sec. IIVI Throughout this paper, we use 
the geometrical units of c = G = 1. 



II. BASIC EQUATIONS 

I — I 

We start by considering the Newtonian limit of a black hole disk system [261 ]. The equations of motion for a test 
particle in this background are very simple. We use the cylindrical coordinates (tzj, tp, z). 

A point mass with a mass M is located at the origin, while a disk exists on the equatorial plane (z = 0) . A smooth 
distribution of disk matter is assumed. If the radial gradient of the density is much smaller than vertical one, we can 
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approximate the density as p 



p(z). A minimal but realistic model for a rotating thick disk may be described by 



Emden's equation [3l| . As in 27[, ignoring the radial gradient, we find that Emden's equation for disk matter density 
p(z) is given by 



^ + *7(7 - 2)p^ 3 (gV = -top, (2.1) 



where k and 7 = 1 + l/n are the polytropic constant and the polytropic index, respectively. The matter density p 
should obey the Poisson equation V 2 Vd = top, where Vd is the potential of the disk. For the isothermal case (7 = 1), 
Eq. (|2.ip has the analytic solution, 

p(z) = -^-sech 2 (-) , (2.2) 



toz V z o 

which corresponds to the disk potential 

Vd(z) — az lncosh [ — ) , (2.3) 
\zoJ 

where zq and a describe the "thickness" of a disk and the surface mass density, respectively. These two parameters 
determine the polytropic constant by the relation 2k = uzq. In the limit of zq — > 0, we recover the potential of an 
infinitcsimally thin disk (Vd ~ a|z|). The corresponding matter distribution is given by the S function from Eq. (|2.2|) . 

Thus, the dynamics of a test particle with a mass p moving around a system of a point mass with a smooth thick 
isothermal disk will be governed by the following (effective) Hamiltonian 

L 2 M , ( z 



H = p 



azo In cosh ( — 



(2.4) 



2 2 2^ 2 tu 2 Vn7 2 + z 2 ' "~ "\z 0/ 
where L is the angular momentum of a particle and the dot denotes the time derivatives. 

III. NUMERICAL ANALYSIS 
A. Two phases of chaos in particle motion 

At first, we analyze particle motion. We numerically integrate the equations of motion of a test particle. The 
symplectic scheme is used because we have the analytic form of the Hamiltonian in this system. The integrated 
time is enough long such that a particle moves thousands times around the central mass. The numerical accuracy is 
monitored by the conservation of the Hamiltonian, which is typically 10 -8 ~ 10 -9 . It guarantees that our numerical 
calculation is reliable. We set M = 1 to fix our units. There are two parameters of a disk which we can change, i.e., 
the surface density a and the width zq. Which parameter dependence we should analyze ? When we change a, there 
are two extreme limits, i.e., a — > and a — > 00, in which the system becomes integrable (see Eq. (|2.4jl ). The gravity 
by the central mass becomes dominant when a — > 0, while the force driven by the disk is dominant as a — > 00. On 
the other hand, if we consider two extreme limits of zq, i.e., the limits of zq — > and of zq — > 00, we find that the 
system is still nonintegrable even in such limits 



2q . Our main aim is to make a distinction between various types of 
chaotic behaviours. For our purpose, the comparison of the cases with different values of zq may not be appropriate. 
Hence we analyze the cases with different values of a, which may provide us continuous change from a regular orbit 
to a very strongly chaotic one. 
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FIG. 1: Poincare maps of orbits of a particle with H = —0.2 and L — 1 across the plane z — in a system of a point mass 
with a disk. We set the thickness the of the disk zo = 0.5, and its surface density (a) a = 0.01, (b) a = 0.1, or (c) a = 10.0. 
All figures are superpositions of trajectories starting from different initial conditions. In Figs, (a) and (c), all trajectories form 
regular tori. In Fig. (b), some trajectories from certain initial conditions still seem to form tori, but others do not. In fact, one 
initial condition, {w, , z, v z ) = (1.2,0,0,0.76) gives an almost two dimensional map on which the orbital points are widely 
scattered, which means that the particle motion is chaotic (We call it Orbit (b)). The outermost trajectories in Fig. (a) and 
(c) are called Orbit (a) and (c), whose initial conditions are (ro, z, v z ) — (1.2,0,0,0.76) and (2.5,0,0,0.49), respectively. 



The parameters of particle orbits such as the energy and angular momentum are appropriately chosen such that 
the motion is bounded. We choose the orbital parameters as H = —0.2 and L = 1, and the disk width as z$ — 0.5. 
Figure [1] shows a set of Poincare maps for different values of the surface density ((a) a = 0.01, (b) a = 0.1, and 
(c) a = 10.0). The equatorial plane (z — 0) is chosen for a Poincare section. We plot the points on the (w,v^) 
plane when the particle crosses the Poincare section with v z > 0. In these figures, trajectories starting from various 
initial conditions are shown. From Figs. [1] (a) and (c), we confirm these system are almost integrable. The outermost 
trajectories are called Orbit (a) and Orbit (c). On the other hand, a widespread chaotic sea is found in Fig. [1] (b). 
This is because the forces by the point mass and by the disk are comparable and those are competing each other. In 
Fig. [T](b), we see the "outermost" trajectory with the initial condition of (tu, z, v z ) = (1.2, 0, 0, 0.76) (called Orbit 
(b)) is not a simple torus but forms an almost two dimensional distribution in which the orbital points are widely 
scattered. It means that the particle motion is chaotic. 

Figure [5] shows the time evolution of the Lyapunov exponents for Orbits (a), (b), and (c) in Fig. [TJ Here we show 
a "local" Lyapunov exponent defined in Appendix [XJ We only refer to the integration time interval t& to define 
it (see Appendix in more details). t& is chosen to be t& — 10 4 , which satisfies the condition of irj -C t& <C tr 
with £_d(« 10 ~ 10 2 ) and £t(~ 10 6 ) being the dynamical time and the total integration period of our calculation, 
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FIG. 2: The time evolutions of "local" Lyapunov exponents for Orbits (a), (b), and (c) in Fig.[T] The Lyapunov exponents for 
Orbits (a) and (c) settle down to very small values, but that for Orbit (b) is large and changing in time. It decreases to a very 
small value in the time interval of t/M = (1.6 ~ 3.8) x 10 5 . 



respectively. We also calculate it with other time intervals, £a = 2 or 4 x 10 4 . We find that the result is not sensitive 
to this choice. We numerically calculate the exponents with the algorithm shown in 29( and show the maximum 
component of it. 

The value for Orbit (a) is very small and almost constant [jjjj]. For Orbit (c), the system is not exactly integrable. 
The system in the limit of a — > oo is of course integrable, but there is no bound orbit in such a limit. Since we are 
analyzing a bound orbit, even if a is very large, we cannot ignore the gravitational effect of a point mass. It makes 
the system nonintegrable. Nevertheless, the motion looks very regular (see the Poincare map in Fig. [T]). In fact we 
find a very small Lyapunov exponent, which is smaller than that of Orbit (a) as shown in Fig. [5] This value is also 
almost constant, which means that the strength of the chaos does not change much in time. Hence we may regard 
this orbit as a regular one. 

On the other hand Orbit (b) gives large positive Lyapunov exponent. It also shows time variation. We should note 
that the value quickly goes down to a very small one in the time interval of t/M = (1.6 ~ 3.8) x 10 5 . We pick up 
the data around this interval and show the time evolution of the r-position of the particle and the Poincare map in 
Fig. [3] From this, we find that although the particle motion in Orbit (b) is chaotic, it stays around r ~ 1.2 — 2.2M 
in the time interval of t/M = (1.6 ~ 3.8) x 10 5 . The motion in this period seems to be nearly regular. In fact, the 
"local" Lyapunov exponent decreases to 5 ~ 6 x 10~ 3 , which is almost the same as those of Orbits (a) and (c). We 
call this phase of motion Orbit (b-2). The phase before this interval, in which a particle motion looks very chaotic, is 
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FIG. 3: The particle motion in the r-direction in terms of time and its Poincare map for the time interval between t/M = 10 5 
and t = 4x 10 for Orbit (b). The Poincare map shows there exist many small tori around the origin. 




FIG. 4: The particle motion of Orbit (b) in the r-direction for a longer time interval than that in Fig. [3] There exists new 
stagnant motion as the same one in Fig. [3] 



called Orbit (b-1). We have also performed numerical integration for longer time period and confirm such phases as 
Orbit (b-2) often appears in a chaotic orbit (see Fig. [4|. The important point is that two different phases of motion 
appear and both a nearly integrable and a more strongly chaotic motion co exist in the same trajectory. 

The Poincare map of Orbit (b-2) in Fig. [3] shows that many small tori exist. It is well known that such a structure 
appears if an orbit is nearly integrable and produces the so-called 1/ f fluctuation [3, [it], 



20 



2l|. Then we also 

analyze the power spectrum of the r-component of Orbit (b-2), which clearly shows a 1/ / fluctuation for / < 10~ 2 M -1 
(see Fig. [5]). This confirms our previous result 2l| in the present model. 



B. Indication of chaos in gravitational waves 



Next we study how to extract information from such a chaotic system and distinguish the orbits, i.e., a nearly 
integrable and more strongly chaotic motions. 

In [21I ] , the authors focused on the power spectrum of particle motion moving in Schwarzschild spacetime and found 
that it shows a power-law behaviour. In this work, we use a similar analysis for the gravitational waves emitted from 
our system. It could be a new and robust way to observe chaotic behaviors in astrophysical objects, as mentioned 
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FIG. 5: The power spectrum of the motion in the r-direction of Orbit (b-2), which is part of nearly regular motion in Orbit 
(b) ( t/M = (1.6 ~ 3.8) x 10°). We find a power-law spectrum, e.g., a 1/ / fluctuation. This is a reflects the existence of small 
tori in the phase space, and the particle moves almost regularly there [21J. 



in our Introduction. The gravitational waves from the present system are calculated by the quadrupole formula 
which is given by 

3 2 6> + l) ,o sin 2 6» 



[{h%-hQ y )cas2ip + 2h% sm2^] 
- (h® z cos ip + h® z sin ip) sin 9 cos 9 , 
[2h% y cos 2cp - 



-/i«)sin2 ¥ >] 



sm if 



h® z cos if) sin ( 



(3.1) 
(3.2) 



where 



h w - = 
i j — 



2d 2 Q l 



dt 2 



with Qij = fj, ( Z % Z 1 - -Sij Z' 



(the reduced quadrupole moment of a point mass). (3.3) 



(r,0,(p) [or (x,y,z)] is the position of a distant observer in spherical coordinates [or Cartesian coordinates], and.Z(t) 
is a trajectory of a particle. We assume that the observer is on the equatorial plane, i.e. (9, if) — (n/2,0). Figured 
shows the waveforms from Orbits (a), (b), and (c). The left panels show the "+" polarization modes of those waves, 
while the right ones are the "x" polarization. The top, middle, and bottom panels correspond to the waves from 
Orbits (a), (b), and (c), respectively. The waves from Orbits (a) and (c) show a periodic feature, which is expected 
from the Poincare maps in Fig.[T] On the other hand, the waves from Orbit (b) show a completely different behaviour. 
We find much random spiky noise in the waveform before t/M = 1.6 x 10 5 and after t/M — 3.8 x 10 5 . This is a typical 



feature of the gravitational waves from highly chaotic motion [11|, |l7J. We also find that the amplitude decreases 
for the time interval of t/M — (1.6 ~ 3.8) x 10 5 . As shown in Fig. [3l in this time interval, the particle moves near 
the small tori in the phase space. This adjective feature of this particle motion appears clearly in the gravitational 
amplitudes. That is, in the phase of a nearly regular motion, the particle position and its velocity do not change 
much compared with those in the more strongly chaotic phase (b-1) (see Fig. [Tfb) and Fig. [3fb)). The time variation 
of the quadrupole moment of the system is small and hence the wave amplitude decreases as well. 

We also calculate the energy spectra of the gravitational waves, which will be one of the most important observable 
quantities in the near future. In Fig. [71 we show the energy spectra for each orbit. Figures [TJa) and (c) show many 
sharp peaks at certain characteristic frequencies. If a motion is regular, we expect several typical frequencies with those 
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harmonics. So such a result reflects that the particle moves regularly. Figure [7] (b) gives the spectrum of Orbit (b). It 
is clearly different from the previous two almost regular cases. It looks just like white noise, below a typical frequency 
f M ~ 10~ 2 , i.e., the shape of the spectrum is flat and it contains many noisy components. However, the spectrum 
of Orbit (b-2) (Fig. EIh-2)), which is analyzed by the orbit only in the time interval of t/M = (1.6 ~ 3.8) x 10 5 , does 
not do so. Rather it looks similar to the spectrum of a regular orbit. Contrary to Fig. [7jb), it does not contain much 
noise at the low frequency region (/M < 10~ 2 ). 

To see more detail, dividing the time interval of Orbit (b) into two, we show the magnifications of the spectra of 
Orbits (a), (b-1), (b-2), and (c) in Fig. [5] Compared to the spectra (a) and (c), the spectra (b-1) and (b-2) contain 
many noisy spikes. Such noisy spikes are usually found in the gravitational waves from a chaotic orbit 17]. However, 
the spectra (b-1) and (b-2) are completely different. The spectrum (b-1) is just white noise. No structure is found. 
On the other hand, the spectrum (b-2) looks similar to those for regular orbits. The "sharp" peaks appear at some 
frequencies, but the widths of those peaks are broadened by many noisy spikes. Therefore, we conclude that Orbit 
(b-2) looks nearly regular but still holds its chaotic character, and such a feature imprints in the spectrum of the 
waves. The important point is that two phases in the particle orbit (b), i.e., the nearly regular phase and the more 
strongly chaotic one, are also distinguishable in the gravitational wave forms and the energy spectra. With this 
analysis, we could constrain orbital parameters. 

IV. SUMMARY AND DISCUSSION 

In this paper we have investigated chaos characteristic for a test particle motion in a system of a point mass with 
a massive disk in Newtonian gravity. To distinguish such characteristics, we propose the gravitational waves emitted 
from this system. At first, we analyzed the motion of the particle by use of the Poincare map and the "local" Lyapunov 
exponent. We found that the phase in which particle motion becomes nearly regular always appears even though the 
global motion is chaotic. We emphasize that both phases of nearly regular and more strongly chaotic motions are 
found in the same orbit. 

The gravitational wave forms and their energy spectra have been evaluated by use of the quadrupole formula in 
each case. In two almost regular cases, the waves show the periodic behaviour and certain sharp peaks appear in those 
energy spectra. In the chaotic case, we have found that the waves show two phases, the nearly regular phase and 
a more strongly chaotic one. In the nearly regular phase, wave amplitude gets smaller in the more strongly chaotic 
phase. The energy spectra are also clearly different. The spectrum in the more strongly chaotic phase looks like white 
noise, but in the nearly regular one, it becomes similar to those in the regular ones. However it is accompanied by 
many small noisy spikes, which is a characteristic feature of a chaotic system. These spikes make the widths of the 
spectrum peaks broader than those in the regular cases. Comparing information from the waves with the particle 
motion, we conclude that we can extract chaotic characteristics of a particle motion of the gravitational waves of the 
system. In the present analysis, in the spectrum (b-2) of the gravitational waves, we do not find a power-law structure, 
which appears in the spectrum of the particle motion. This may be because the waveform is given by the change of 
the quadrupole moment, which contains higher time derivatives of a particle trajectory such as acceleration. It may 
be much more interesting if one can find the 1// behaviour in some information of the gravitational waves because 
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FIG. 6: The gravitational waveforms evaluated by the quadrupole formula. Top, middle, and bottom figures correspond to 
those for Orbits (a), (b), and (c), respectively. The left and right rows give the "+" and "x" polarization modes, respectively. 



such an indication may specify the type of chaos more clearly. This is under investigation. 

Finally, we mention a possibility to constrain parameters in a dynamical system. If the gravitational waves are 
observed for a sufficiently long time, we can monitor the time variation of the wave amplitudes, their forms and 
polarizations. We can then calculate the energy spectra for some durations. If the spectra show one of the typical 
characteristics found in this paper, the parameters of a particle motion could be constrained. Of course, a realistic 
system can be more complicated, and the present model may be too simple. But we believe the characteristic behaviour 
of the gravitational waves found in this paper will help us to understand a chaotic system. Therefore our next task 



is to anal yze the gravitational waves from various chaotic systems, especially relativistic chaotic systems 

3, 0, mill Q, H E 
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32[. Then, we should investigate whether or not the correlation between the gravitational 
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FIG. 7: The energy spectra of the gravitational waves shown in Fig. [6] Orbit (b-2) gives the spectrum of the waves for the 
"stagnant motion", i.e., when the particle motion of Orbit (b) becomes near regular for t/M = (1.6 ~ 3.8) x 10 J . Figures (a) 
and (c) show many sharp peaks at certain characteristic frequencies. This is because of the regular motion. The spectrum in 
Fig. (b), which looks like white noise for fM < 10 -2 , is clearly different from those in Figs, (a) and (c), but the spectrum 
in Fig. (b-2) does not look like white noise. It looks similar to the cases (a) and (c). However, the peaks are not sharp 
but rather broadened by appearing so many other spikes. Note that the typical frequency of the orbits is in the range of 
fM = 10" 2 ~ 10" 1 (see Fig. ©. 



waves and chaos in dynamical systems found in this work is generic. 
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APPENDIX A: LOCAL LYAPUNOV EXPONENT 

In this appendix, we give the definition of "local" Lyapunov exponent. Our definition of "local" Lyapunov exponent 
is somewhat different from the conventional one[3|, but those are essentially the same. 

At first, let us consider the system whose time evolution is described by a set of differential equations in N- 
dimensional space, 

x = F(x), (Al) 

where x(i) is a iV-dimcnsional vector. 

The time evolution of the orbital deviation (5x, which is the difference between two nearby orbits, obeys the following 
set of linear differential equations: 

*x=^(x(t))Jx. (A2) 

The solution of Eq. (|A2|) can be written formally as 

<5x(t) = E/* to <5xo, (A3) 

where <5x is an "initial" deviation at some time to and U l tQ is an evolution matrix, which is given by the following 
integration; 



U\ = exp 



(A4) 



We define the "local" Lyapunov exponent in time interval [to,t] by 



, 1 \\U\ ei A U\ e 2 A • • • A U\ eJ 

A(e , t) = log " t0 ,, *2-= to_*H A5 

t-to ||ei Ae 2 ---Ae fe || v y 

for fc = 1, 2, • • • , N, where e k is a /c-dimensional subspace in the tangent space at the initial point Xo, which is spanned 
by k independent vectors ej (i — 1,2, •• • ,k), A is an exterior product, and || o || is a norm with respect to some 
appropriate Riemannian metric. If we take a limit of t — > oo, \(e k ,oo) correspond to the conventional Lyapunov 
exponents. 

If the integration time interval t& = t — to is much longer than the dynamical time of the system, we may find 
convergent values for each X(e k ,t), which are almost independent of t a (or to)- We may call them "local" Lyapunov 
exponents at t. The maximum value of "local" Lyapunov exponents, i.e. A(f) = max{A(e fc , t)\k = 1, 2, • • • , N} is the 
most important one for our discussion. So we also call it the "local" Lyapunov exponent at t. 



8 




1 1 1 1 1 1 1 1 

-6 -5 -4 -3 -2 -1 
Log 10 [fM] 



